## "Where are the missing dead" replication
##
## Reproduce: Monthly outcome plot
##			  - Figure 2 (16 month outcome around the announcement)
##
## 			  
## 
##
## Input:     prov_month.dta 
##			  - 16 months of data, 4291 observations


## Guoer Liu
## 7/30/2022 


setwd("~/Dropbox/Missing dead/replication/")


library(haven)
library(extrafont) ## change fonts in plot



### Function 

dataCollapse <- function(vector, year.v, treat.ind){
	# Args:
	#	vector: 	n x 1 vector to be collapsed by province-treat
	#	year.v: 	n x 1 vector
	#   treat.ind:  n x 1 vector of treatment indicator
	out <- tapply(vector, 
				  list(year.v, treat.ind), 
				  function(x) mean(x, na.rm = T))
	return(out)
}


paraPlot <- function(collapMat, 
					matrix.x.lab = NULL, matrix.title = NULL,
					l.mar = 3.5,
					hline = NULL,
					x.axis.vec = NULL,
					y.lab = TRUE, y.lab.names = NULL,
					y.min = NULL, y.max = NULL,
					legend = TRUE, leg.pos = NULL){
	# Args:
	#	collapMat:      n x 2 matrix (column 1: control; column 2: treat)
	#   matrix.x.lab:   a string of figure x axis label
	#	matrix.title:	a string of figure title
	#	hline:			a scalar indicate the position of vertical line
	#	x.axis.vec:		a length n string vector that shows x axis
	#   y.lab:			logical parameter to turn off y label
	#   y.lab.names:    a string of y axis lable
	#   y.min:			a scalar indicates the minimum value on y axis
	#   y.max:			a scalar indicates the maximum value on y axis
	#   legend:			logical parameter to turn off legend
	#	leg.pos:		a string indicates legend position
	pch.cex <- 1.6
	ctl.v <- collapMat[, 1]
	tr.v <- collapMat[, 2]
	#y.max <- max(x)
	#y.min <- min(x)
	x.len <- nrow(collapMat) 

	par(mar= c(6, l.mar, 3, 1))
	plot(NA, xlim = c(0.5, x.len), 
	         ylim = c(y.min, y.max), 
	         type = "n", 
	         xlab = matrix.x.lab,
	         main = matrix.title, 
	         ylab = "", axes = F, cex.main = 1.5, cex.lab = 1.3)
	grid(nx = NULL, ny = NULL,
     	 lty = 2,      # Grid line type
     	 col = "gray", # Grid line color
     	 lwd = 0.6)
	abline(v = hline, col = alpha("orangered", 0.8), lwd = 1.5)
	lines(ctl.v, lty = 2, lwd = 1.5)
	points(ctl.v, cex = pch.cex, pch = 1, col = "black")
	lines(tr.v, lty = 1, lwd = 1.5)
	points(tr.v, cex = pch.cex, pch = 19, col = "black", bg = "black")
	if (y.min < 1){
		text(seq_len(x.len), 0, labels = x.axis.vec, srt = 20, adj = 1, 
			xpd = TRUE)
	} else {
		text(seq_len(x.len), 1.33, labels = x.axis.vec, srt = 20, adj = 1, 
			xpd = TRUE)
	}
	axis(1, at = seq(1, x.len, 1), label = rep("", x.len), tick = TRUE)
	#	 cex.axis = 1, las = 1)
	if (y.lab == TRUE){
		axis(2, at = seq(ceiling(y.min), floor(y.max), 1), 
			label = seq(ceiling(y.min), floor(y.max), 1),
			cex.axis = 1, las = 1)
		title(ylab = y.lab.names, line = 2, cex.lab = 1)
	}
	if (legend == TRUE){
		legend(leg.pos, title = NULL, 
			   legend = c("Treated", "Control"), 
			   col = "black", 
			   lty = c(1, 2), cex = 0.8,
			   pch = c(19, 1),
			   lwd = c(1, 1))
	}
}







mydata.m <- read_dta("prov_month.dta")


deathsum_pmlg_tr <- dataCollapse(mydata.m$deathsum_pmlg_tr, 
								 mydata.m$month, mydata.m$treat_group)
totcount_pmlg_tr <- dataCollapse(mydata.m$totcount_pmlg_tr, 
								 mydata.m$month, mydata.m$treat_group)
traf_deathsum_pmlg_tr <- dataCollapse(mydata.m$traf_deathsum_pmlg_tr, 
									  mydata.m$month, mydata.m$treat_group)
traf_totcount_pmlg_tr <- dataCollapse(mydata.m$traf_totcount_pmlg_tr, 
									  mydata.m$month, mydata.m$treat_group)
othr_deathsum_pmlg_tr <- dataCollapse(mydata.m$othr_deathsum_pmlg_tr, 
									  mydata.m$month, mydata.m$treat_group)
othr_totcount_pmlg_tr <- dataCollapse(mydata.m$othr_totcount_pmlg_tr, 
									  mydata.m$month, mydata.m$treat_group)

mon.names <- c(NA, "2003-08", NA, "2003-10", NA, "2003-12", NA, "2004-02", NA, 
			  "2004-04", NA, "2004-06", NA, "2004-08", NA, "2004-10") 

#pdf("f2_monthOutcome.pdf", family = "Times New Roman",
#	 width = 10, height = 6.18)
par(mfrow = c(2, 3), mgp=c(3.5, 0.8, 0))
paraPlot(deathsum_pmlg_tr, "Time", "Aggregate",
		hline = 8, 
		x.axis.vec = mon.names, 
		y.lab = TRUE, y.lab.names = "Log (Deaths)",
		y.min = 1.8, y.max = 5.0,
		legend = FALSE, leg.pos = "bottomright")
paraPlot(traf_deathsum_pmlg_tr, "Time", "Traffic",
		hline = 8, 
		x.axis.vec = mon.names, 
		y.lab = FALSE, y.lab.names = "Log", l.mar = 1.5,
		y.min = 1.8, y.max = 5.0,
		legend = FALSE, leg.pos = "bottomright")
paraPlot(othr_deathsum_pmlg_tr, "Time", "Others",
		hline = 8, 
		x.axis.vec = mon.names, 
		y.lab = FALSE, y.lab.names = "Log", l.mar = 1.5, 
		y.min = 1.8, y.max = 5.0,
		legend = TRUE, leg.pos = "topright")
paraPlot(totcount_pmlg_tr, "Time", "Aggregate",
		hline = 8, 
		x.axis.vec = mon.names, 
		y.lab = TRUE, y.lab.names = "Log (Acc. counts)",
		y.min = 0.5, y.max = 4.0,
		legend = FALSE, leg.pos = "bottomright")
paraPlot(traf_totcount_pmlg_tr, "Time", "Traffic",
		hline = 8, 
		x.axis.vec = mon.names, 
		y.lab = FALSE, y.lab.names = "Log", l.mar = 1.5,
		y.min = 0.5, y.max = 4.0,
		legend = FALSE, leg.pos = "bottomright")
paraPlot(othr_totcount_pmlg_tr, "Time", "Others",
		hline = 8, 
		x.axis.vec = mon.names, 
		y.lab = FALSE, y.lab.names = "Log", l.mar = 1.5, 
		y.min = 0.5, y.max = 4.0,
		legend = TRUE, leg.pos = "topright")
#dev.off()